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Abstract 

The calculation of the solvation properties of a single water molecule in liq- 
uid water is carried out in two ways. In the first, the water molecule is placed in 
a cavity and the solvent is treated as a dielectric continuum. This model is an- 
alyzed by numerically solving the Poisson equation using the DelPhi program. 
The resulting solvation properties depend sensitively on the shape and size of 
the cavity. In the second method, the solvent and solute molecules are treated 
explicitly in molecular dynamics simulations using Ewald boundary conditions. 
We find a 2 kcal/mole difference in solvation free energies predicted by these 
two methods when standard cavity radii are used. In addition, dielectric con- 
tinuum theory assumes that the solvent reacts solely by realigning its electric 
moments linearly with the strength of the solute's electric field; the results of 
the molecular simulation show important non-linear effects. Non-linear solvent 
effects are generally of two types: dielectric saturation, due to solvent-solute 
hydrogen bonds, and electrostriction, a decrease in the solute cavity due to 
an increased electrostatic interaction. We find very good agreement between 
the two methods if the radii defining the solute cavity used in the continuum 
theory is decreased with the solute charges, indicating that electrostriction is 
the primary non-linear effect and suggesting a procedure for improvement of 
continuum methods. The two methods cannot be made to agree when the 
atomic radii are made charge independent, but charge dependent cavity radii 
are shown to greatly improve agreement. 
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Introduction 

Dielectric continuum methods based on the numerical solution of the Poisson 
Equation provide estimates of solvation free energies. Because this approach is much 
faster than full molecular dynamics simulations, continuum theory (CT) has been 
applied to a wide variety of systems, from aqueous solutions of large solute molecules 
such as proteins and nucleic acids to small solutes such as atomic ions [?]. The 
widespread use of such continuum calculations makes it important to test them 
against simulations with explicit solvent molecules and against experimental data 
where available. Previous comparisons of continuum theory and molecular simula- 
tions have found good agreement for solvation energies, typically about 10% [?, ?]. 
The focus of this note will be to compare continuum theory with full molecular simu- 
lations of a small polar solute, for equilibrium properties, such as solvation energies, 
electrostatic potentials, and average electric fields, from which the electrostatic forces 
can be calculated. 

The solvation process studied here is the solvation of a single water molecule 
in liquid water. The free energy of solvation, AA so i, of a water molecule can be 
decomposed into a three-step process: 1) removal of the gas-phase partial charges from 
the gas-phase water molecule (AAf° s ); 2) insertion of the uncharged water molecule 
into the solvent making a hydrophobic cavity in the bulk water (AA cav ); and 3) 
charging the water solute to the desired liquid phase charges (AA es ) (see Figure |T|) [?]. 
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The free energy, AAf" s , for the first step is the self-energy for electronic polarization 
in the gas-phase or the free energy difference between the isolated molecule with its 
gas phase partial charges and with no partial charges. The free energy for the liquid- 
state charging step, AA es , contains two parts, a self-energy for polarization in the 
liquid phase and a solvent electrostatic contribution. The two polarization terms will 
not completely cancel since the gas phase and solution charges are not necessarily (or 
generally) equal and also the polarization energy will depend on the medium. The 
two polarization energies are commonly neglected [but see [?, ?]]. 

The solvent electrostatic contribution to the free energy is the focus of this paper 
and is the subject of comparison between the full molecular simulation and the con- 
tinuum calculations. This step consists of reversibly charging the hydrogen charge, 
Qh, on the solute water molecule from to 0.5e (the oxygen charge, Qo, equals - 
2Q H ) while keeping the solvent charges equal to the usual liquid state values. (We 
are including charges higher than Q^ 9 — higher than is necessary for the evaluation 
of AA es for H 2 — in order to look at the effects of high charge.) In the continuum 
approximation, all the solvent electrostatic properties are contained in the dielectric 
constant, e. Setting e equal to the same value given by the molecular potential used 
in the simulation, we can assure self-consistency. 

The properties of interest are the energy of the solute, the average electrostatic 
potential and the field of the solvent at the solute charge sites. The electrostatic 
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potential, 0, at position ri Q (where 1 labels the solute molecule and a labels the 
charge site) due to the solvent is given by 

N 3 
3=2 13=1 

The electrostatic potential energy of the solute is 

(U„> = Qa(<Pa) (2) 

a=l 

and the electric field is 

-<W~(^) (3) 
where the brackets (■ • •) indicate an average over the solvent configurations. The free 
energy, AA es , for the solute with a set of partial charges Q can be calculated using 
the charging integral, 
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AA es (Q) = J2Q a / 1 (0(r lQ )) A dA 







(XJ es ) x /\d\ (4) 



where the angular brackets, (■ ■ -)\, indicate an ensemble average with the solute 
charges equal to AQ [?]. Eq.([|) is integrated numerically by performing simulations 
at 11 different hydrogen charges, ranging from to .5. 

Continuum theories assume that the solvent response to the solute is linear in the 
solute charge. Specifically, the potential, 0(r lQ ), is linear in the charge Q a and there- 
fore, from Eq.(f|), AA es oc Q 2 . The results of the simulations will test the assumption 



of linear response. Specifically, continuum theory based on the Poisson equation is a 
linear theory in that it assumes that the orientations of the solvent molecules respond 
linearly to the electric field of the solute [?]. In real molecular liquids, non-linear 
responses not described by the Poisson equation arise from many factors including 
the formation of solute-solvent hydrogen bonds and electrostriction, the decrease in 
the excluded volume of the solute due to the increased solute-solvent Coulombic in- 
teraction. Strong hydrogen bonds prevent the solvent from further responding to 
the solute causing dielectric saturation which decreases the solvation energy. On the 
other hand, electrostriction increases the solvation energy. 

Methods 

The water potential used here is the SPC potential, characterized by an OH bond 
length of 1 A, an HOH bond angle of 109.47°, charges on the hydrogens and oxygen 
equal to 0.41e and -0.82e, respectively, and a Lennard- Jones interaction between 
oxygen atoms with a well depth, e/ks, equal to 78.2 Kelvin and a radius, a, equal to 
3.166 A[?]. The molecular dynamics simulations were performed on the Connection 
Machine CM-5 using a direct N 2 method[?], with 511 solvent molecules and 1 solute 
molecule. Periodic boundary conditions, using the Ewald sum for the long-ranged 
electrostatic potentials, a time step of 1 femtosecond and the SHAKE algorithm for 
enforcing bond constraints are used[?]. The simulations are done in the canonical 
(constant T,V,N) ensemble by coupling to a Nose thermostat [?, ?] and are at a 



5 



density of 1 g/cm 3 and a temperature of 300 Kelvin. Each data point at a given 
charge represents a 40 picosecond simulation [?]. The electrostatic potential with 
periodic boundary conditions is 

N 

4>(ria) = E / EE c Wl r i« _i V/3 + n l ( 5 ) 

n j=l p 

where the prime on the sum over periodic images n indicates that for n = the term 
j=l is omitted. In the standard Ewald evaluation of Eq. (|5]), the energy is written as 
a sum of four terms, 

N 

0(r lQ ) = Q/3erfc(/t \r la - r jP \)/ \r la - r jf3 \ 

i=2 p 

N A 1 

+ E 7T 2 e- G2/4K2 cos(G • (r lQ - r,,)) 

j=l (3 *-* G^O U 



N 



31/3 3=1 P 

~ EQ/3 erf ( K l r lQ - r ip\)/ \ r la ~ r ip\ (6) 

P 

where k is parameter in the Ewald sum chosen for computational convenience to be 
6/L, L is the length of the primary simulation cell, and G is a recriprocal lattice 
vector of the periodic simulation cells [?]. 

In the continuum calculations, the solute is characterized by a molecular cavity 
defined by spheres around each charge site and a dielectric continuum outside this 
cavity. The potential, 0, can then be found by solving the Poisson equation 

V- e(r)V0(r) +47rp(r) = (7) 



where e(r) is a position dependent dielectric constant, equal to 1 inside the solute 
cavity and 65 outside (65 is the dielectric constant of the SPC water model [?, ?]), 
and p(r) is the charge density. Equation ([7]) is solved using the DelPhi program, 
which discretizes space on a cubic grid (with 65x65x65 points) [?]. This approach is 
thus based on a finite-difference solution of the Poisson Equation. The calculations 
reported here use focusing boundary conditions, in which successive calculations are 
done, each with a finer mesh, using the previous coarser grid results to correct for 
the long range boundary effects. At the highest resolution, 80% of the grid points are 
inside the solute cavity, corresponding to a grid spacing of about 0.034 A. The water 
molecule geometry is the SPC geometry. The oxygen cavity radius, To, is set equal 
1.77 A and the hydrogen cavity radius is 0.8 A. The oxygen radius is about 2 -5/,6 cr, 
one half the minimum of the Lennard- Jones potential. This is the standard method 
for choosing radii in DelPhi calculations [?, ?]. 

Results and Discussion 

The free energies corresponding to the molecular dynamics (md) simulations and 
to the CT calculations are shown in Figure |2|. At the charge value of SPC water 
(Qi/=-41), the md yields -8.4±.5 and CT yields -10.5 kcal/mole for AA es . The re- 
ported error bars are two standard deviations. This difference of about ± 2 kcal/mole 
between the two methods is comparable to the agreement between the CT and molec- 
ular free energy calculations (using the TIP4P water potential) for several solutes (but 
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not water) [?]. The charge dependence for the simulation data is not quadratic, as 
can be seen by the poorness of a quadratic fit (the dotted line). 

Previous calculations for SPC and also TIP4P water embedded in a dielectric con- 
tinuum find AA es =-10.96 and -10.89 kcal/mole, respectively [?]. These calculations 
use an oxygen cavity radius of 1.5A and a hydrogen cavity radius of 1.16A. The sensi- 
tivity of the CT results to these input parameters will be discussed below. Other CT 
calculations by Sharp, et al. [?], for the TIP4P geometry, which include polarizability 
of the solute, find AA es =-9.3 kcal/mole. There have been calculations of AA so i from 
molecular simulations, but the electrostatic part, AA es , was not reported so there is 
no other molecular simulation to compare to [?]. 

Figure ^ shows the average electrostatic potentials, {(f)). The md and CT results 
are qualitatively different in two respects: 1) {(f>n) and {4>o) from the md simula- 
tions are not linear in charge and 2) {(f) h) and {4>o) from md are not zero at zero 
charge, whereas the 0's from the CT calculations are linear in charge and zero at 
zero charge. As is well known, molecular water solvates a small uncharged solute by 
forming a clathrate cage around the uncharged spheref?]. For realistic water poten- 
tials, the electrostatic potential is not zero inside the cage because fluctuations of 
the water molecules that form the clathrate cage are known to bring the hydrogen 
atoms closer to the cage than the oxygen atoms. This is apparent from the charge 
distribution function (see below). Since the positive charges can get closer than the 
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negative charges there is a net electrostatic potential inside the sphere. However, it 
is the spatial dependence and not the exact value of <fi that is important because 1) 
properties such as the electric field and the energy are invariant with respect to an 
additive constant in cf) and 2) the exact value of ((f)) is strongly dependent on the 
choice of boundary conditions used in the Ewald sum. In the Ewald summation, 
the periodically replicated system must be surrounded by a medium, in Eq. (FJ) the 
medium is taken to be an insulator (this boundary term is the third term on the 
right-hand side of the equation). If the surrounding medium is taken to a conductor, 
then this term vanishes[?]. The boundary conditions have only a slight influence on 
most properties, however, for the electrostatic potential the boundary conditions give 
rise to a large constant term. 

The electric field, E, arising from the solvent at each solute site is shown in Figure 
f|. The geometry of the water molecule is shown in the upper left corner of Figure f|. 
The top panel shows the x-component of the field, E x , at the position of the hydrogen 
Hi, E x at the other hydrogen is minus the field at Hi and E x at the oxygen is zero. 
The bottom panel shows the y-component of the field at the oxygen and hydrogen 
positions. Since there is no net force on the molecule and for the continuum solvent 
all forces are Coulombic then 

£Q Q E(n Q ) = 0, 

a 

the y-component of the field at the oxygen site and hydrogen sites must be equal. 
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For the SPC simulations, there is an additional force on the oxygen atom due to 
the non-bonded Lennard- Jones interaction which has a small nonzero component in 
the y-direction (at Q#=0.41, it is 1.47 kcal/mole/A) so the net force is zero but the 
electrostatic fields do not exactly balance. The field from the md simulation is not 
linear in charge, unlike the CT field. This is consistent with the simulation results 
for the electrostatic potential, <ft, being nonlinear and the free energy, AA es , being 
non-quadratic. (Figures |3|-(|). A non-quadratic charge dependence was reported by 
Jayaram, et al[?], for the solvation of a spherical cation. 

Any conclusions regarding the validity of dielectric continuum theory are of course 
dependent on the input parameters, particularly on the radii, io and r#. The values 
for the AA e <,, the x- and y-components of the electric field and the Coulomb energy, 
(U e5 ), are given in Table [I]. (From the fact that the electrostatic potential (0) is 
linear in charge, it follows from Eqs. |2] and |]that AA es = (U cs ) /2 for the continuum 
results.) Also shown on Table [I] are the SPC simulation results and the results of two 
approximate models (see below). Some cavity radii give accurate estimates for the free 
energy, but it is not possible for continuum theory to simultaneously predict values 
of AA es , (U es ), and the electric field, E, to all within 20%. In addition, the charge 
dependence of (U es ) and E must be linear and AA es must be quadratic, in contrast 
to the results of the simulation. There are sets of radius parameters which will give 
the best quadratic fit to the simulation free energy. Two such sets are (r G =1.88 A, 
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r#=0.80 A) and (ro=l-83 A, r#=1.00 A) which will give the free energy curve shown 
by the dotted line of Figure 

The solvent reorganization that results from increasing the charge of the solute 
molecule is reflected in the radial distribution of charge, defined by 

47T 3 

Q a (r) = — pr 2 £ Q/J&tfto ( 8 ) 

6 13=1 

where p is the bulk density of the solvent, the factor of three in the dominator is 
introduced to reflect that there are three atomic sites, and g a p is the pair correlation 
function between the atomic sites a on the solute and (3 on the solvent molecules [?]. 
Figure |5| shows Qo( r ) — the radial distribution of charge as a function of distance 
from the solute oxygen site — for three values of the solute charge: Q#=0,0.25, and 
0.50. The distribution functions oscillate between positive hydrogen atom peaks and 
negative oxygen atom peaks. As the charges on the solute are increased, the first 
peak moves in. In addition, a peak grows in at a O-H hydrogen bond distance of 
1.8 A. This hydrogen-bond peak is faintly visible at Q#=0.25 and is prominent at 
Q ff =0.50. 

We focus on two of the possible explanations of the breakdown of dielectric con- 
tinuum models: 1) dielectric saturation, in which the orientational ordering (from 
hydrogen bonds) of the first solvation shell decreases the dielectric response of the 
solvent and 2) electrostriction, in which the solvent molecules come into closer contact 
with the solute molecule as its polarity is increased, and which thus leads to a de- 
ll 



crease in the solute cavity as solute charge is increased. Dielectric saturation effects 
decrease the solvation energy whereas electrostrictive effects increase the solvation 
energy. Both of these solvent responses are visible in Figure [5| We will take two 
approaches to understanding the differences in the md simulation free energies and 
the continuum free energies. In one, we will scale the cavity size with solute charge 
in the continuum model. In the other, we will explicitly include first shell waters in 
the CT calculations. 

Scaled Radius. From Figure ^j, for Q#=0 the cavity is larger than 2~ 5 / 6 cr (=1.78 
A), since the solvent peak does not start until after 2 A, and the CT results, which 
use 1.77 as a cavity radius, therefore overestimate the free energy at low charges (see 
Figures 0). This raises a question as to whether the differences between the simulation 
and CT calculations are due to the approximation of a continuum solvent or just an 
inconsistent choice of cavity size. Estimates of the cavity size can perhaps be found 
from an examination of the liquid structure, although assigning a sharp solute/solvent 
boundary from a continuous distribution is ambiguous. A simple heuristic method is 
to find the radius which gives the best value for the Coulomb energy, (U es ). For each 
value of Qh the optimal oxygen radius is found (the hydrogen radius is set equal to a 
constant value of 0.8 A) and the free energy is calculated from Eq. |[ The resulting 
(U es ) and AA es are in almost exact agreement with the simulation values (see Table 
p. Of course, this agreement is by construction and with one adjustable parameter 
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it is a trivial accomplishment to fit the simulation energies. However, this method 
also gives good values for the electric fields, which is not by construction (see Figure 
^[). The optimized oxygen radius is shown for each value of the hydrogen charge in 
Figure ^ (there is no optimized value at Q#=0 since at this point, (U es )=0). The 
oxygen radius shows a strong charge dependence — the radius varies by 30% over the 
range of charges — and this dependence is approximately linear with a slope of about 
1.2 A/e. There is some uncertainty in the radius for low values of the charge since 
the energy is small and variations of ± .1 A in tq give rise to energies which are all 
within the error bars of the simulation. The values of the scaled tq are shown by the 
arrows on Figure^ indicating that these radii are consistent with the liquid structure. 
The arrow on the Q#=0 plot is the radius for Q#=0.05. 

Semicontinuum Methods. In this model the nearest neighbor solvent water molecules 
are treated explicitly and the rest are treated as a continuum. This method could 
include both dielectric saturation effects as the first solvation shell orders around the 
solute and electrostriction effects as the first solvation shell moves closer to the solute. 
The focus of this analysis is to see if we can explain the non-linear effects seen in the 
simulation in terms of reorganization of the local solvation shell. This method is 
implemented as follows. Configurations generated from the SPC simulation are taken 
and the coordinates of the solute plus its n nearest neighbors are inputed into the 
DelPhi program. This is done for about 300 configurations, each taken .1 picoseconds 
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apart. The energy and fields on the central solute will now be a sum of a part coming 
from the n explicit nearest neighbors and the dielectric continuum. This model is 
similar to previous semicontinuum methods, although here 1) we are averaging over 
configurations of the first shell molecules and not using thermodynamic data for the 
first shell contribution as is done in the other studies, and 2) we are taking an ar- 
bitrary surface for regions beyond the first shell (as defined by the radii of the first 
shell atoms) and not a spherical shell [?, ?]. Another semicontinuum study by Rashin 
and Bukatin includes one first shell water inside the continuum and thermodynamic 
properties are found by integrating numerically over the solute and single solute water 
coordinates and then extrapolating to a full hydration shell[?]. We are treating the 
number of explicit water molecules, clS db variable, ranging from 4 to 8, in order to 
measure how many solvent molecules are strongly influenced by the solute, 

Figure |8| compares the results of the full simulation, the continuum theory (with 
ro=1.77 A and r#=0.80 A) and the semicontinuum results with 4 and 8 neighbors. 
This figure gives AA es divided by Qh versus Qh- The semicontinuum results and the 
full simulation exhibit a similar non- linear dependence on AA es /Q# on Qh, although 
the deviations from linearity are not as great for the semicontinuum model. This 
means that the non-linear effect is largely local, as also supported by the scaled radius 
analysis. At low charges the n=4 semicontinuum results are essentially the same as the 
pure continuum results, implying that we have not included enough explicit neighbors. 
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At high charges the semicontinuum results are different from the continuum results. 
The free energy at a given charge depends on values of the potential at lower charges, 
through Eq. (|4j). Continuum theory with radii ro=1.77A and r#=0.8A overestimates 
(Ues) at low charges (because the oxygen radius is too small, see the preceding section) 
and underestimates it at higher charges (see Table [1]). The n=4 semicontinuum 
results also overestimate the potential energy at low charges (because it has not 
included enough neighbors and is similar the pure continuum results) and at higher 
charges, it does not have the same fortuitous and compensating underestimation of 
{U es ) (because at higher charges four solvating waters can more adequately describe 
the local environment). Therefore, the continuum theory gives better estimates of 
the free energy at high charge than the n=4 semicontinuum method. For other 
properties, such as (U es ) (Table [[]) and the electric fields (Table [I] and Figure ||), 
the semicontinuum methods are closer to the simulation results and as n increases, 
the agreement improves. At Q#=.41 the first solvation shell contains 5 neighbors 
as measured by integrating the goo out to the first minimum. At Q#=0, the first 
solvation shell is broader and contains more molecules (about 16). By including up 
to 8 explicit neighbors in the semicontinuum model we are then including all the 
first shell solvent molecules of the higher charge solute, but only some of the first 
shell molecules of the lower charge solute. At the highest charges, above .41, the 
semicontinuum results are the same for 4 and 8 neighbors, so at these charges there 
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is strong tetrahedral ordering 4 nearest neighbors are enough to describe the local 
structure. 

Conclusion 

The simple example of solvation presented here is that of a single water molecule 
with partial charges on the atomic sites which are varied between zero and values 
greater than the charges for water (a partial charge on the hydrogen atoms equal to 
0.50e). The continuum theory (CT) calculation by definition shows a simple quadratic 
dependence on charge for the free energy, AA es , and a linear dependence for electro- 
static potential, <fi and the electric fields, E. From the md simulations, the charge 
dependence of AA es is definitely non-quadratic and that of (0) and E non-linear (see 
Figures 0, |3] and ^). The md simulations then indicate that some type of solvent 
reorganization other than the re-orientation of the electric moments of the solvent 
molecules is occurring as the solute is charged. Because these effects are important 
and because they are not included in the CT calculations, continuum theory cannot 
simultaneously predict free energies, potential energies and electric fields to within 
20% no matter what radius parameters are used (see Table |l]). There are sets of 
radius parameters which will predict the free energies shown by the dotted line in 
Figure [| which represents the best quadratic fit to the simulation free energy. 

Two simple models are suggestive of the types of solvent reorganizations that are 
important to this solvation process. These models also provide a method to improve 
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continuum theory. In one model, the size of the solute cavity is charge dependent 
(scaled radius) and in the other model the system inside the dielectric continuum 
contains the solute together with explicit neighboring solvent molecules (semicontin- 
uum). The energies and free energies for all of the different methods (md simulation, 
pure continuum theory, scaled radius, and semicontinuum) are summarized in Table 
|I[ As can be seen from this Table and also from Figure [], the scaled radius method 
provides a good estimate of the free energy and the electric fields. The agreement 
between the scaled radius and md simulation results for the Coulomb energy, (U es ), is 
by construction since (U es ) is used to determine the scaled radius. The success of the 
scaled radius calculations suggest that electrostriction plays a more important role 
than dielectric saturation of close-lying solvent molecules in explaining the deviations 
of the continuum model from the full molecular solvent. The oxygen radius used in 
the scaled radius method varies by a large amount (30%) with the charges of the 
solvent (see Figure |7]). This large decrease in the solute cavity is also seen in the 
liquid correlation functions from the simulation (see Figure |5[). 

In conclusion, simple continuum theory, as usually implemented in DelPhi calcula- 
tions, is not capable of determining free energies to better than 2 kcal/mole accuracy 
and gives inaccurate predictions of electrostatic potentials and electric fields. We be- 
lieve that continuum theory can be improved to give better than 1 kcal/mole accuracy 
by adopting charge dependent radii to allow for electrostrictive effects. It remains to 



17 



invent a theoretical model which predicts how site radii should be scaled with charge. 
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Table 1: Comparison of the md simulation results (first row) with continuum theory 
(CT) results for a variety of cavity radii (jo, *h), the scaled radius CT results and 
semicontinuum results with n explicit neighbors. The properties listed are the free 
energy, AA es , the Coulomb energy, (U es ), both in kcal/mole and the x-component, 
E x , and y-component, E^, of the electric field in kcal/(mole e A) (see Figure 4) for a 
solute with Q#=0.41. The numbers in parenthesis are two standard deviation error 
estimates. 





AA es 


(u es ) 


Ei(rHi) 




Simulation 


-8.4(5) 


-23.32(6) 


22.2(1) 


33.6(1) 


CT (r o =1.77,r H =0.80) 


-10.5 


-21.0 


25. 


32. 


CT (r o =1.88,r H =0.80) 


-8.2 


-16.4 


17. 


25. 


CT (r o =1.83,r H =1.00) 


-8.2 


-16.4 


14. 


25. 


CT (r o =1.80,r H =1.00) 


-8.6 


-17.2 


15. 


26. 


CT (r o =1.70,r H =1.10) 


-8.7 


-17.4 


12. 


27. 


CT (r o =1.50,r H =1.16) 


-10.1 


-20.2 


12. 


31. 


Scaled radius CT 


-8.4 


-23.6 


26. 


35. 


semicontinuum(n=4) 


-11.3(3) 


-26.1(1) 


26.1(7) 


39.2(4) 


semicontinuum(n=6) 


-10.3(4) 


-25.2(4) 


24.0(7) 


37.6(5) 


semicontinuum(n=8) 


-9.7(5) 


-24.2(2) 


22.9(3) 


36.1(2) 
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Figure 1: Thermodynamic cycle for the solvation of water in water. 

Figure 2: Free Energy, AA es , from md simulations (solid line) and DelPhi continuum 
calculations (dashed line) as a function of solute charge (in kcal/mole). The dotted 
line is a quadratic fit to the simulation data. 

Figure 3: The averags electrostatic Potential, (0), in kcal/mole/e, at the hydrogen 
(solid lines) and oxygen (dashed lines) position, comparing the SPC simulations (the 
lines showing data points and error bars) with DelPhi continuum results (no point 
symbols). 

Figure 4: Electric field, E, in kcal/mole/e/A in the x-direction at the H\ hydro- 
gen atom (top), and in the y-direction at the hydrogen and oxygen atoms (bottom) 
comparing md simulation (solid (H) and dotted lines(O)) and DelPhi continuum cal- 
culations (dashed lines). In the continuum model E y at the H and O sites are equal, 
for the simulation there is a small difference (see text). See the upper left-hand corner 
for the definition of the coordinate system. 

Figure 5: Solvent charge distribution (determined from the md simulations) about the 
solute oxygen site for a) Qh=0, b) Q H =0.25 and c) Q H =0.50. The arrows indicate 
the optimized oxygen cavity radius. 

Figure 6: As for Figure 4, but comparing md simulation (solid lines), DelPhi contin- 
uum results with fixed radii (dashed lines), and scaled radii continuum results (dotted 
lines). 

Figure 7: Optimized oxygen cavity radius, To, as a function of solute charge. 

Figure 8: Free energy divided by hydrogen charge, Qh, for the full md simulation 
(long-dashed line), semicontinuum with 4 neighbors (dotted line), semicontinuum 
with 8 neighbors (short-dashed line), and the continuum theory (solid line). 

Figure 9: As for Figure 4, but comparing simulation (solid lines) results with semi- 
continuum results for 4 (dotted lines) and 8 neighbors (dashed line). 
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